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Abstract 

The purpose of this tutorial is to introduce the main concepts behind normal and 
anomalous diffusion. Starting from simple, but well known experiments, a series 
of mathematical modeling tools are introduced, and the relation between them is 
made clear. First, we show how Brownian motion can be understood in terms of a 
simple random walk model. Normal diffusion is then treated (i) through formalizing 
the random walk model and deriving a classical diffusion equation, (ii) by using 
Pick's law that leads again to the same diffusion equation, and (iii) by using a 
stochastic differential equation for the particle dynamics (the Langevin equation), 
which allows to determine the mean square displacement of particles, (iv) We discuss 
normal diffusion from the point of view of probability theory, applying the Central 
Limit Theorem to the random walk problem, and (v) we introduce the more general 
Fokker-Planck equation for diffusion that includes also advection. We turn then 
to anomalous diffusion, discussing first its formal characteristics, and proceeding 
to Continuous Time Random Walk (CTRW) as a model for anomalous diffusion. 
It is shown how CTRW can be treated formally, the importance of probability 
distributions of the Levy type is explained, and we discuss the relation of CTRW 
to fractional diffusion equations and show how the latter can be derived from the 
CTRW equations. Last, we demonstrate how a general diffusion equation can be 
derived for Hamiltonian systems, and we conclude this tutorial with a few recent 
applications of the above theories in laboratory and astrophysical plasmas. 
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1 Introduction 



The art of doing research in physics usually starts with the observation of a 
natural phenomenon. Then follows a qualitative idea on " How the phenomenon 
can be interpreted" , and one proceeds with the construction of a model equa- 
tion or a simulation, with the aim that it resembles very well the observed 
phenomenon. This progression from natural phenomena to models and math- 
ematical prototypes and then back to many similar natural phenomena, is the 
methodological beauty of our research in physics. 

Diffusion belongs to this class of phenomena. All started from the observations 
of several scientists on the irregular motion of dust, coal or pollen inside the air 
or a fluid. The roman Lucretius in his poem on the Nature of Things (60 BC) 
described with amazing details the motion of dust in the air, Jan Ingenhousz 
described the irregular motion of coal dust on the surface of alcohol in 1785, 
but Brownian motion is regarded as the discovery of the botanist Robert 
Brown in 1827, who observed pollen grains executing a jittery motion in a fluid. 
Brown initially thought that the pollen particles were "alive", but repeating 
the experiment with dust confirmed the idea that the jittery motion of the 
pollen grains was due to the irregular motion of the fiuid particles. 

The mathematics behind "Brownian motion" was first described by Thiele 
(1880), and then by Louis Bachelier in 1900 in his PhD thesis on "the theory 
of speculation" , in which he presented a stochastic analysis of the stock and 
option market. Albert Einstein's independent research in 1905 brought to the 
attention of the physicists the main mathematical concepts behind Brownian 
motion and indirectly confirmed the existence of molecules and atoms (at that 
time the atomic nature of matter was still a controversial idea). As we will see 
below, the mathematical prototype behind Brownian motion became a very 
useful tool for the analysis of many natural phenomena. 

Several articles and experiments followed Einstein's and Marian Smoluchowski's 
work and confirmed that the molecules of water move randomly, therefore a 
small particle suspended in the fluid experiences a random number of impacts 
of random strength and direction in any short time. So, after Brown's obser- 
vations of the irregular motion of "pollen grains executing a jittery motion", 
and the idea of how to interpret it as "the random motion of particles sus- 
pended inside the fluid", the next step is to put all this together in a flrm 
mathematical model, "the continuous time stochastic process." The end re- 
sult is a convenient prototype for many phenomena, and today's research on 
" Brownian motion" is used widely for the interpretation of many phenomena. 

This tutorial is organized as follows: In Sec. 2, we give an introduction to 
Brownian motion and classical random walk. Sec. 3 presents different models 
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for classical diffusion, the Langevin equation, the approach through Pick's law, 
Einstein's approach, the Fokker-Planck equation, and the central limit theo- 
rem. In Sec. 4, the characteristics of anomalous diffusion are described, and a 
typical example, the rotating annulus, is presented. Sec. 5 introduces Contin- 
uous Time Random Walk, the waiting and the velocity model are explained, 
methods to solve the equations are discussed, and also the Levy distributions 
are introduced. In Sec. 6, it is shown how, starting from random walk models, 
fractional diffusion equations can be constructed. In Sec. [7] we show how a 
quasi-linear diffusion equation can be derived for Hamiltonian systems. Sec. 8 
briefly comments on alternative ways to deal with anomalous diffusion. Sec. 
9 contains applications to physics and astrophysics, and Sec. 10 presents the 
conclusions. 



2 Brownian Motion and Random walks 

2. 1 Brownian Motion Interpreted as a Classical Random Walk 

To build a firm base for the stochastic processes involved in Brownian motion, 
we may start with a very simple example. 
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Fig. 1. Random walk in one dimension (along the vertical axis) as a function of time 
(to the right on the horizontal axis). 

We consider a random walk in one dimension (ID) and assume that the parti- 
cles' steps are random and equally likely to either side, left or right, and of 
constant length I (see Fig. [T]). The position z^r of a particle starting at = 
after N steps is 

N 

zn = Az^ + Azn_i + -I- Azi = Azi, (1) 

i=l 
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so that the squared length of the path equals 



^2 . 



' N \ / ^ \ ^ 

,?=i / \fc=i / j,fc=i 



N N N 

= Az|+ ^ZjAzk = Nf+ J2 ^zjAzk. (2) 

When averaging over a large number of particles, we find the mean squared 
path length as 

<zl>=Ne + l Az,Az,\. (3) 

\i,fc=i,jVfc / 



Each step of the walk is equally likely to the left or to the right, so that 
the displacements Azi are random variables with zero mean. The products 
AzjAzk are also random variables, and, since we assume that Azj and Az^ 
are independent of each other, the mean value of the products is zero, so that 
the expectation value of the mixed term in Eq. ([3]) is zero. We thus find 

< >= Ne. (4) 



The root-mean square displacement after steps of constant length I (mean 
free path) is 



R:=J<zl> = iVN. (5) 



We can now estimate the number of steps a photon starting from the Sun's 
core needs to reach the surface of the Sun. From Eq. ([5]), we have N = {R/iY 
and since the Sun's radius is ~ 10^'' cm and the characteristic step (taken 
into account the density in the solar interior) is ~ 1 cm, we conclude that 
photons make 10^° steps before exiting from the Sun's surface (this can answer 
questions like: if the Sun's core stops producing energy, how long will it take 
until we feel the difference on Earth?). 

The mean free path i can be estimated with a simple model. By assuming 
that a particle is moving inside a gas with a mean speed < v >, the distance 
traveled between two successive collisions is ^ =< f > r, where r is called 
collision time. If the particle has radius a and travels a distance L inside the 
gas with density n, then it will suffer Aira'^Ln collisions, which is just the 
number of particles in the volume 47ra^L the particle sweeps through. The 
mean free path is then defined through the relation Aira'^in = 1, i.e. i is the 
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distance to travel and to make just one collision, so that 



We may thus conclude that the number of steps a particle executes inside a 
gas during a time t is N = t/r, and, with Eq. (j4]) and the above relation 
i =< V > T, the mean squared distances it travels is 

< >= Nf = {t/T){< v>T)e = {<v> i)t. (7) 

Assuming that the random walk takes place in 3 dimensions and that the gas is 
in equilibrium and isotropic, we expect that < >=< y"^ >=< >= ^ ^ , 
and the mean square path length in 3 dimensions is 

< >= 3 < V > it = Dt, (8) 

where Z^:=3<t;>^is called the diffusion coefficient, which is a useful 
parameter to characterize particle diffusion in the normal case (see Sec. |3]). 
Important here is to note the linear scaling relation between < > and time 
t. 



2.2 Formal Description of the Classical Random Walk 



More formally, we can define the classical random walk problem as follows. 
We consider the position r of a particle in 1, 2, or 3 dimensional space, and 
we assume that the position changes in repeated random steps Ar. The time 
At elapsing between two subsequent steps is assumed to be constant, time 
plays thus a dummy role, it actually is a simple counter. The position fn of a 
particle after n steps, corresponding to time t„ = nAt, is 

fn = Ar; + Afn-l + Afn-2 + •.• + Afi + Tq (9) 



where tq is the initial position, and Afi is the ith step (or increment, or dis- 
placement). The position r„ as well as the increments Af are all random 
variables. To specify the problem completely, we have to prescribe the proba- 
bility distribution q^ff{Af) for the increments Af, which yields the probability 
for the particle to make a certain step Ar, (with given length and direction). 
Writing q^ff{Af) in this way, we have made the assumptions that all the in- 
crements have the same probability distribution and that the increments are 
independent of each other (the value Af takes in any realization is completely 
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independent of the value taken in the previous step by Afi^i). Generahza- 
tions to time-dependent increment distributions or correlated increments are 
of course possible. 

Since r„ is a random variable, the solution we are looking for in the random 
walk problem is in the form of the probability distribution P(r, t„), which 
yields the probability for a particle to be at position r at time t = tn = nAt. 

If we are interested in the mean square displacement, we can just square Eq. 
(Q, and, rearranging the terms in the same way as in Eq. ([2]), we find for 
ro = 

n n 

{r'n)= E E i^r-An). (10) 

j=l,{k=j) j,k=l,kj^j 



The first term on the right hand side is just a sum over the variances cr^^j 
of q/\^{Ar), since by definition cr^^ := {Ar^) = J Af'^q^f^^Af) dAr if the 
mean value of q^if{Ar) is zero. The second term on the right hand side is the 
covariance cov(Ar^-, Afk) of the random walk steps, and it is zero if the steps 
a particle takes are independent of each other. We thus can write Eq. ffTOl) as 

n n 

{r^)= E ^lr,j+ E cov(Ar-;-,Af,). (11) 

j=l,{k=j) j,k=l,k^j 

The particular random walk we considered in Sect. I2.11 can thus be understood 
on the base of Eq. ([TT]) as a case with zero covariance and variance c^r*, j — 
due to the constant step length, which leads to the mean square displacement 
in Eq. dH). 



3 Models for Normal Diffusion 



3.1 Langevin's Equation 



We turn now to a different way of treating Brownian motion. We again consider 
a particle with mass m performing a random walk inside a fluid due to the 
bombardment by the fluid mol ecules, which obe y an equilibrium (Maxwellian) 
distribution. Pierre Langevin (ILangevinl . Il908l ) described this motion with a 
simple but very interesting stochastic differential equation (let us work in one 
dimension for simplicity). 



mx 



-ax + F{t), 



(12) 
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where the term ax represents the friction force, x is the particle velocity, a is 
the damping rate and depends on the radius of the particle and the viscosity of 
the fluid, and F{t) is a random fluctuating force due to the random bombard- 
ment of the particle by the fluid molecules. If the random fluctuating force 
were absent, the particle starting with an initial velocity vq would gradually 
slow down due to the friction term. Multiplying Eq. f|T2l) with x, we have 



mxx 



m 



d{xx) 
dt 



— X 



—axx + xF{t), 



and after taking averages over a large number of particles we find, since < 
xF{t) >= due to the irregular nature of the force F{t), 



d < XX > 



m- 



dt 



m < X > —a < XX > . 



(13) 



Since the background gas is is in equilibrium, the kinetic energy of the particle 
is proportional to the gas temperature, m < x^ > /2 = kT/2, where k is the 
Boltzmann constant and T the temperature of the gas. Eq. (1131) now takes the 
form 



d 



kT 



— + 7 <xx>= — , 
,dt I m 



where 7 = a/m, which has the solution 



< XX >= ; = Ce ^ H . 

2 dt a 



(14) 



At t = 0, the mean square displacement is zero, so that 
Eq. (IT^ becomes 



C + kT/a, and 



1 d < x2 > kT 



dt 



— (1 
a 



On integrating the above equation we find the solution 
2kT 



< X > -- 



'--(1 

7 ^ 



— e 



-7t 



(15) 



In the limit t << I/7 (time much shorter than the collision time) the solution 
in Eq. (fT5|) is of the form < >~ (expanding the exponential up to 
second order), which is called "ballistic" diffusion and means that at small 
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times particles are not hindered by collisions yet and diffuse very fast, see 
Sect. 14.21 In the other limit, t » I/7, the solution has the form 

< >~ 1, (16) 



or, for the 3 dimensional case, if again the gas is in equilibrium and isotropic 
so that < > /3 =< >, 

< >= 1 = Dt, (17) 



where D = 6kT/a is an expression for the diffusion constant in terms of 
particle and fluid characteristics (cf. Eq. (IHl)), and note that again < > has 
a simple scaling relation with time, < >= Dt, as in Sect. 12.11 



3.2 Modeling Diffusion with Pick's Law 



Diffusion usually occurs if there is a spatial difference in concentration of 
e.g. particles or heat etc., and it usually acts such as to reduce the spatial 
inhomogeneities in concentration. 

Let us consider particle diffusion along the ^-direction in 3-dimensional space, 
and let us assume that two elementary areas perpendicular to the flow (in 
the x-y-plane) are a distance Az apart. Particle conservation implies that the 
time variation of the density n{z, t) inside the elementary volume AxAyAz 
equals the inflow minus the outflow of particles, so that, if J{z, t) denotes the 
particle flux, 

Oti( z t] dJ 

— AxAyAz = J{z)AxAy - J{z + Az)AxAy = -—AxAyAz 



which leads to the diffusion equation in its general form, 

dn{z,t) dJ{z,t) 

= d^- ^^^^ 

The problem that remains is to determine the particle flux J. From its physical 
meaning, it obviously holds that 

J{z,t) = n{z,t)v{z,t), (19) 
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where v{z,t) is an average particle flow velocity. Using this expression in Eq. 
(1181) leads to a closure problem, we would need to flnd ways to determine 
v{z,t). 

It is well documented experimentally that the flux of particles J crossing a 
certain area (again, say in the x-y-plain) is proportional to the density gradient 
along the z-axis (Pick's Law), 

J, = -Diz)^, (20) 



where D is the diffusion coefficient discussed already in the previous sections, 
and which generally may also depend on z. With Eq. (120!) . the diffusion equa- 
tion takes the classical form 

dniz^ = #D(.)M!ll) , (21) 
dt dz ^ ' dz ' ^ ' 



or, for constant diffusion coefficient, 

dn{z,t) _ ^ d^n{z,t) 

dt ~ dz^ • ^ ^ 



In infinite space, and if all particles start initially from z = 0, the solution of 
Eq. ([22D is 



n{z,t) = -^e-'/'^\ (23) 



where Nq is the total number of particles inside the volume under considera- 
tion. The solution obviously is identical to a Gaussian distribution with mean 
zero and variance 2Dt. The variance is defined as 

< z^(t) >= I z^n{z, t) dz = 2Dt, (24) 



which is just identical to the mean square displacement, so that the results 
obtained earlier, using the simple version of the random walk in Sect. l2.1l or the 
stochastic differential equation of Langevin in Sect. 13. H are again confirmed. 

Diffusion obeying Eq. (12^ is called normal diffusion and is characteristic 
for the diffusion processes in systems that are in equilibrium or very close to 
equilibrium. Generalizing the above results (Eqs. (l2T]) . (l22l) . ( l23l) . (l24l) 
is simple and can be found in the literature (see references). 
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3.3 Einstein's Formalism for the Classical Random Walk and the Diffusion 
Equation 



A different approach to treat norm al diffusion was introduced by Bachelier 
and by Einstein (see 



Einstein! . Il905l ). Here, tlie starting point is tlie classical 
and we consider the 1-dimensional case. 



random walk as defined in Sec. 12.21 
According to Sec. 12.21 the solution of the random walk problem is in the form 
of the probability distribution P{z, t) for a particle at time t to be at position 
z. Assume that we would know the distribution P{z,t — At) one time-step 
At earlier (remember At is assumed constant). If particles are conserved, the 
relation 



P{z, t) = P{z - Az,t- At) qAz{Az) 



(25) 



must hold, with q^^ the distribution of random walk steps. Eq. (!25|) states that 
the probability to be at time t at position z equals the probability to have 
been at position z — Az at time t — At, and to have made a step of length 
Az. We still have to sum over all possible Az, which leads to the Einstein (or 
Bachelier) diffusion equation, 

oo 

P{z, t)= j P{z - Az,t- At) qAz{Az) dAz (26) 

oo 



This is an integral equation that determines the solution P{z, t) of the random 
walk problem as defined in Sec. 12. 2[ The power of this equation will become 
clear below when we will show ways to treat cases of anomalous diffusion. Here, 
we still focus on normal diffusion. As will become clear later, it is actually a 
characteristic of normal diffusion that the particles take only small steps Az 
compared to the system size. This implies that qAz{Az) is non-zero only for 
small Az, the integral in Eq. fl26|) is only over a small Az-range, and we can 
expand P{z — Az, t — At) in z and t (also At is small), 

P{z -Az,t- At) = P{z, t) - At dtP{z, t) - Az dzP{z, t) 

+ \Az'dlP{z,t) (27) 

Inserting into Eq. we find 



P{z,t)= j P{z,t)qAz{Az) dAz- j AtdtP{z,t)qAz{Az) dAz 
- j AzdzP{z, t) qAz{Az) dAz 
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+ 1 ^ Az^ dlP{z, t) q^z ( Az) dAz (28) 

Obviously, P{z,t), its derivatives, and At are not affected by the integration, 
furthermore we use the normahzation of qaz (/ QAzi^z) dAz = 1), we assume 
qaz to be symmetric (mean value zero, / Az q^ziAz) dAz = 0), and we use 
the definition of the variance (/ Az"^ q^ziAz) dAz = <j\^), so that 



Piz, t) = Piz, t) - At dtP{z, t) + d',P{z, t) 



(29) 



or 



d,Pz,t=^^d!P{z,t) (30) 



i.e. we again recover the simple diffusion equation, as in Sec. 13. 2[ with diffu- 

2 

sion coefficient D = and the solution to it in infinite space is again the 



2At ' 

Gaussian of Eq. 



3.4 Fokker- Planck Equation 



The Fokker-Planck (FP) equation (or Kolmogorov forward equation) is a more 
general diffusion equation than the simple equations introduced in Sees. I3.2l an 
13.31 We again start from a description of diffusion in terms of a random walk, 
as in Sec. 13. 3^ we though relax two assumptions made there: (i) We assume 
now that the mean value /xa^ of the random walk steps can be different from 
zero, which corresponds to a systematic motion of the particles in the direction 
of the sign of /iaz, and (ii) we assume that both the mean and the variance 
can be spatially dependent, /zaz = fJ'Az{z) and cta^ = o''Azi^)y which means 
that the distribution of increments depends on the spatial location, i.e. it is 
of the form qAz,z{Az, z). To be compatible with these assumptions, Eq. (12^ 
must be rewritten in a slightly more general form, 

00 

P(z, t)= J Piz - Az,t- At)gAz,z(A2, z - Az) dAz, (31) 



which is the Chapman- Kolmogorov equation, and where now q/\z,z{Az, z) is 
the probability density for being at position z and making a step Az in time 
At. The FP equation can be derived in a way similar to the one presented in 
Sec. 13.31 We expand the integrand of Eq. (ETj) in a Taylor-series in terms of z, 
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so that P{z, t) = AB dAz, with 

A = P{z, t) - dtP{z, t)At - d,P{z, t)Az + ]^dlP{z, t)Az^ + (32) 



where we have also expanded to first order in t, and which is of course the 
same as Eq. fl27|) . and newly we have 

B = qAzA^^^ ^) - d.q^zA^z, z)Az + ]^dlq^,A^z, z)Az^ + ... (33) 

(note that the Taylor expansion is with respect to the second argument of 
Q'Az.z, we expand only with respect to z, not though with respect to Az). In 
multiplying and evaluating the integrals, we use the normalization of q^z,z 
(/ qAz,z{^^i ^) dAz) = 1), the definition of the mean value 
ifJ'Aziz) := / Azq/^.^^z{^z, z) dAz) and of the second moment {{Az^){z) : = 
J Az^q^z,z{^2:, z) dAz), and expressions like / Az dzqAz,z{.^z, z) dAz are con- 
sidered to equal dz J AzqAz,z{Az, z) dAz = dzfiAz{z), so that, on keeping all 
terms up to second order in Az, we find the Fokker-Planck equation, 

dtP{z,t) = -dz[V{z)P{z,t)] + dl[D{z)P{z,t)], (34) 



with V{z) = fiAziz)/At a drift velocity, and D{z) = ( Az^)(z)/2At the diffu- 



sion coefficient (for a 3 dimensional formulation see e.g. iGardinerl . |2004| ). The 
basic difference between the FP equation and the simple diffusion equation in 
Eq. (130|) is the appearance of a drift term, and that both the drift velocity 
and the diffusion coefficient are allowed to be spatially dependent (Pick's law 
also allows a spatially dependent diffusion coefficient, see Eq. (pT!) ). These 
differences allow the PP equation to model more complex diffusive behaviour. 

The PP equation is also applied to velocity space, e.g. in plasma physics in 
order to treat collisional effects, or to position and velocity space together. It 
has the advantage of being a deterministic differential equations that allows 
to describe the evolution of stochastic systems, as long as the diffusivities and 
drift velocities are known, and as long as the conditions for its applicability 
are met, see the remarks below. 

We can illustrate the typical use of the PP equation on the example of Brown- 
ian motion in the Langevin formalism in Sec. 13.11 which allowed us to calculate 
the diffusion coefficient in Eq. flTTl) . If we are interested in the evolution of the 
distribution of particles P, then with the Langevin formalism we would have 
to follow a large number of individual particles over the times of interest and 
then to construct the distribution function, which may become very intense 
in computing effort. Instead, one can use the diffusivity from Eq. (fT7|) . insert 
it into the PP equation, and solve the PP equation. Since in this example the 
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diffusion coefficient is constant and there is no drift velocity, the FP equation 
reduces to the simple diffusion equation (122|) . whose solution for P is given in 
Eq. 



We just note that in the general case where the diffusion coefficient is z- 
dependent, D = D{z), there is a small difference between the diffusive term 
in the Fokker-Planck equation and the Fickian diffusion equation, Eq. fl2ip . 
in that the diffusivity is only once di fferentiated in the latter. This difference 



and its consequences are analyzed in Ivan Milligen et al.l (120051 ) 



/^From its derivation it is clear that the FP equation is suited only for sys- 
tems close to equilibrium, with just small deviations of some particles from 
equilibrium, or, in the random walk sense, with just small steps of the parti- 
cles performing the random walk, exactly as it holds for the simple diffusion 
equation in Sect. 13.31 

A further natural generalization for a diffusion equation in the approach fol- 
lowed here would be not to stop the Taylor expansion in Eqs. (1521) and (1551) 
at second order in z, but to keep all terms, which would lead to the so-called 

Kramers- Moyal expansion. 

Mor e details about t he Fokker-Planck equation can be found in the literature 



e.g. iGardinerl . |2004J ) . 



3.5 Why Normal Diffusion Should Be the Usual Case 



The appearance of normal diffusion in many natural phenomena close to equi- 
librium and the particular Gaussian form of the solution of the diffusion equa- 
tion can also be understood from probability theory. The Central Limit Theo- 
rem (CLT) states that if a statistical quantity (random variable) is the sum of 
many other statistical quantities, such as the position of a random walker after 
n steps according to Eq. (Q, and if (i) all the Azi have finite mean /xaz and 
variance cr\^, (ii) all the Azi are mutually independent, and (iii) the number 
n of the Az is large, then, independent of the distribution of the A^j's, the 
distribution P{z,tn) of Gaussian. In particular, if /i^.^ = 0, 2:0 = and 

all the Azi have the same variance, then 

1 - 

P(z,t„) = - -e ^ (35) 



with variance a^^ = na\^, or, if we set n = t/ At, then = t„cr^^/At. 
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The mean square displacement equals per definition the variance, 

(zitnY) = J z^P{z, tn) dz = = tncrU At (36) 

(for zq = 0), and diffusion is thus always normal in the cases where the CLT 
applies. 

Moreover, the CLT predicts quantities that are the result of many small scale 
interactions to be distributed according to a Gaussian, and indeed this is 
what we found for the distribution of the classical random walker, see Eq. 
0231) . Stated differently, we may say that the appearance of non-Gaussian 
distributions is something unexpected and unusual according to the CLT. 
We just mention that also the equilibrium velocity distributions of gas or 
fluid particles are in accordance with the CLT, the velocity components, say 
Vx, Vy, Vz, follow Gaussian distributions, and therewith the magnitude v = 
+ Vy + v1 exhibits a Maxwellian distribution. Again then, the appearance 
of non-Maxwellian velocity distributions is unexpected on the base of the CLT. 



4 Anomalous Diffusion 

4-1 Systems Far from Equilibrium: The rotating Annulus 



video camera 




1 1 o 



d d 2 d 

Fig. 2. Rotating annulus. 
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The simple experiment of the rotating annulus, shown in Fig. p allows to illus- 

trate the differences between norrnal and anomalous diffusion ( ISolomon. Weeks fc Swinney 
I994J : IWeeks. Urbach fc Swinneyl . Il996l ). Water is pumped into the annulus 
through a ring of holes marked with / and pumped out through a second 
ring of holes marked with O. The annulus is completely filled with water and 
rotates as a rigid body (the inner and outer walls rotate together). The pump- 
ing of the fluid generates a turbulent flow in the annulus. A camera on top of 
the annulus records the formation of the turbulent eddies inside the rotating 
annulus and allows to track seeds of different tracer particles injected into the 
fluid and to monitor their orbits (see Fig. [3]). 




Fig. 3. (a) The formation of eddies inside the rotating annulus, as recorded by the 
camera (left panel), and (b) typical orbits of tracer particles inside the annulus 
(right panel). 



In the case of normal diffusion, which occurs mainly in fluids close to equilib- 
rium, the particle trajectories are characterized by irregular, but small steps, 
which makes trajectories look irregular but still homogeneous (see Fig.[T]). The 
trajectories shown in Fig. [3] for the highly turbulent rotating annulus, which 
is far away from equilibrium, show different types of orbits, with two basic 
new characteristic, there is "trapping" of particles inside the eddies, where 
particles stay for 'unusually' long times in a relatively small spatial area, and 
there are "long flights" of particles, where particles are carried in one step 
over large distances, in some cases almost through the entire system. 



4-2 The Scaling of "Anomalous" Trajectories 



Normal diffusion has as basic characteristic the linear scaling of the mean 
square displacement of the particles with time, (r^) ~ Dt. Many different 
experiments though, including the one shown in the previous section, reveal 
deviations from normal diffusion, in that diffusion is either faster or slower, 
and which is termed anomalous diffusion. A useful characterization of the 
diffusion process is again through the scaling of the mean square displacement 
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with time, where though now we are looking for a more general scaling of the 
form 



(At)) ~ t\ 



(37) 



Diffusion is then classified through the scaling index 7. The case 7 = 1 is nor- 
mal diffusion, all other cases are termed anomalous. The cases 7 > 1 form the 
family of super-diffusive processes, including the particular case 7 = 2, which 
is called ballistic diffusion, and the cases 7 < 1 are the sub-diffusive processes. 
If the trajectories of a sufficient number of particles inside a system are known, 
then plotting log < > vs log t is an experimental way to determine the type 
of diffusion occurring in a given system. 

As an illustration, let us consider a particle that is moving with constant 
velocity v and undergoes no collisions and experiences no friction forces. It then 
obviously holds that r = vt, so that {r^(t)) ~ t^. Free particles are thus super- 
diffusive in the terminology used here, which is also the origin of the name 
ballistic for the case 7 = 2. Accelerated particles would even diffuse faster. 
The difference between normal and a anomalous diffusion is also illustrated in 
Fig. m where in the case of anomalous diffusion long "flights" are followed by 
efficient "trapping" of particles in localized spatial regions, in contrast to the 
more homogeneous picture of normal diffusion. 



Fig. 4. (a) Random walk in dynamical systems close to equilibrium (normal diffusion; 
trajectory on the left), (b) random walk in dynamical systems far from equilibrium 
(anomalous diffusion; trajectory on the right). 

It is to note that anomalous diffusion manifests itself not only in the scal- 
ing of Eq. (137|) with 7 7^ 1 (which experimentally may also be difficult to be 
measured), but also in 'strange' and 'anomalous' phenomena such as 'uphill' 
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diffusion, where particles or heat diffuse in the direction of higher concentra- 
tion, or the appearance of non-Maxwellian distributed particle velocities (see 
Sec. 13. 5p . very often of power-law shape, which is very common in high energy 
astrophysics (e.g. cosmic rays), etc. 



5 Continuous Time Random Walk 

5. 1 Definition 

Given the experimental ubiquity of anomalous diffusion phenomena, the ques- 
tion arises of how to model such phenomena. One way of tackling it is through 
the random walk formalism. So-far, we have used the random walk to model 
classical diffusion, and in Sec. 13.31 it had been shown how the random walk is 
related to a simple diffusion equation if the steps the particles take on their 
walk are small. One way to model anomalous diffusion is by relaxing the latter 
condition, and to allow the particles to also take large steps, where 'large' in a 
finite system means large up to system size, and in infinite systems it means 
that the steps are unbounded in length. Useful in this context is the family 
of Levy distributions as step-size distributions gA^- They are defined in closed 
form in Fourier space (see Sec. 15.3.31 below), and they have the property that 

g^f (A^) ~ |A^|-^-", for |Az| large, < a < 2, (38) 

so that there is always a small, though finite probability for any arbitrarily 
large steps size. The Levy distributions all have an infinite variance, ^ = 
/ Az^g^'^(Az) c/Az = oo, which makes their direct use as a step size distribu- 
tion in the classical random walk of Sec. 12.21 and Eq. ([9]) impossible: Consider 
the case of a random walk in 1-D, with the position of the random walker after 
n steps given by the 1-D version of Eq. ([9]), and the mean square displacement 
(for zq = 0) given by Eq. (ITT]) . Let us assume that the steps are independent of 
each other, so that the covariances are zero and the mean square displacement 
is (z^) = nal^, which is infinite, already after the first step. 

A way out of the problem is to release time from its dummy role and make 
it a variable that evolves dynamically, as the walker's position does. In this 
way, infinite steps in space can be accompanied by an infinite time for the step 
to be completed, and the variance of the random walk, i.e. its mean square 
displacement, remains finite. The extension of the random walk to include the 
timing is called Continuous Time Random Walk (CTRW). Its formal definition 
consists again of Eq. ([9]), as described in Sec. 12.21 and, moreover, the time at 
which the nth step of the walk takes place is now also random (a random 
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variable), and it evolves according to 

tn = Atn + Atn-1 + Atn-2 + •.• + Ati + to, (39) 



where to is the initial time, and the Ati are random temporal increments. To 
complete the definition of the CTRW, we need also to give the probability 
distribution of the Ati, i-e. we must specify the probability for the ith step to 
last a time Ati. 

Two case are usually considered (not least to keep the technical problems 
at a manageable level), (i) In the waiting model, the steps in position and 
time are independent, and one specifies two probabilities, one for Af, the gAr 
already introduced, and one for At, say q^t- Here then At is interpreted as 
a waiting time, the particle waits at its current position until the time At 
is elapsed, and then it performs a spat ial step Af during which no time is 
consumed (e.g. MontroU &: Weiss . 19651 ). (ii) In the velocity model, the time 



At is interpreted as the traveling time of the particle. At = | Ar|/f , where v is 
an assumed constant velocity (the velocity dynamics is not included, usually, 
see though Sec. 15.41) . so that the distribution of increme nts is qAz,At — ^(At — 



Af\/v)qAriAr) (e.g. IShlesinger. West fc Klafted . 119871 ). We just note that in 



the general case one would have to specify the joint probability distribution 
9A2,At(A-2, At) for the spatial and temporal increments. 



5.2 The CTRW Equations 



The CTRW equations can be understood as a generalization of the Einstein 
equation, Eq. (126!) . or the Chapman-Kolmogorov equation, Eq. (1311) . It is useful 
to introduce the concept of the turning-points, which are the points at which a 
particle arrives at and starts a new random walk step. The evolution equation 
of the distribution of turning points Q{z, t) (here in 1-D) follows basically from 
particle conservation. 



Q{z,t)= J dAz j dAtQ{z - Az,t - At)qAz,AtiAz, At) 



+ S{t)P{z,t = 0) + S{z,t), (40) 

where the first term on the right side describes a completed random walk 
step, including stepping in space and in time, the second term takes the initial 
con dition P(z,t = 0) in t o acco unt, and the third term is a source term (see 



e.g. IZumofen fc Klafteii . Il993f ). 



The expression for P{z, t), the probability for the walker to be at position z at 
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time t, is different for tlie waiting and for the velocity model, respectively. In 
case of the waiting model, where gA2,At(A-2, At) = gAt(At)gA2(A-2;), we have 

t 

Pwiz,t) = J dAtQ{z,t- At)^wiAt), (41) 





with ^w(At) := f^ l dt'q AtW) the probability to wait at least a time At (e.g. 
Zumofen fc Klafteii . Il993i l 



In the velocity model, where gAz,At(Az, At) = 6{At—\Az\/v) qAz{Az), P{z,t) 
takes the form 



vt t 



Py{z,t)= J dAz J dAtQ{z- Az,t- At),^viAz,At) 



(42) 



-vt 



with 



oo oo 



$y(A2,At) = h{\Az\-vAt) I dz' j dt'5{t' -\z'\/v)qAz{z' 



(43) 



lA^I At 



the pro bability to make a step of length at least \Az\ and of duration at least 
At (e.g. IZumofen fc Klaftei Il993l : IShlesinger. West fc Klafte^] . Il987l ). 



Both, the expression for Py/ and Py determine the probability for seeing the 
particle when moving in-between two turning points, taking into account only 
the part of the random walk in which time is consumed by the particle. 



The kind of diffusion that the CTRW formalism yields depends on the dis- 
tribution of step increments. If the increments are small, then the treatment 
of Sec. 13.31 can be applied again, diffusion is normal, and again a simple dif- 
fusion equation can be derived. If the increments are not small, then super- 
as well as sub-diffusion can result, depending on the concrete choice of incre- 
ment distributions. For instance, small spatial steps in combination with Levy- 
distributed, long waiting times will yield sub-diffusion in the waiting model. 
An important property of the CTRW equations is that they are non-local, 
both in space and time (which is also termed non-Markovian). Anomalous 
diffusion phenomena in the CTRW approach are thus considered non-local 
processes, and with that they are far from equilibrium processes. 
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5.3 Treating the CTRW Equations 



5.3.1 Remarks on the solution of the CTRW equations 



A standard way to treat the CTRW equation is by transforming them to 
Fourier (F) and Laplace (L) space, whereby the convolution theorems of the 
two transforms are used. We will illustrate this procedure on the example of 
the waiting model below in Sec. 15.3.41 

The CTRW equations are though not always of a convolution type, e.g. the 
velocity model has not a convolution structure anymore, due to the appearance 
of time in the integration limits (see Eq. (H2l) ). so that Fourier Laplace methods 
are not directly applicable anymore, and other methods are needed (actually 
also in the expression for Q, Eq. pUj) . time appears in the Az-integration 
limits in case of the velocity model). 

FL transforms, if applicable, usually do not allow to calculate the probability 

P{z,t) in closed analytical form, but rather some asym ptotic properties of it, 

such as the mean square displacement at l arge times (e.g. lKlafter. Blumen &: Shlesinger 
19871 : iBlumen. Zumofen fc Klafterl . Il989l ). On the other hand, FL transforms, 
if applicable, allow to transform the CTRW equations into other kinds of equa- 
tions , e.g. in one instead of the two i ntegr al equations (IKlafter. Blumen fc Shlesingerl . 
19871 : iBlumen. Zumofen fc Klafteii.ll989l. e.g.). into an integro-differe ntial equa- 



19871 ) 



or even 



tion (or master equation, e.g. IKlafter. Blumen fc Shlesingerl . 
into a fractional diffusion equation, which has the form of a diffusion equations, 
the fractional derivatives are though generalized, non-local differentiation op- 
erators. In Sec. 16.21 we will show how a fractional diffus ion equation arises 
naturally in the context of the waiting model (see also e.g. iMetzler fc Klafter . 

2nnnl . l2nn4h . 



Another standard wa y of treating the CTRW e quati ons is with Monte Carlo 
simulations (see e.g. IVlahos. Isliker fc Lepreti l2004l ). o r else, the eq uations 
can be solved numerically, with an appropriate method (llslikerl . l2008l ). 



5.3.2 Fourier and Laplace transforms of probability densities 

For any probability density function (pdf) such as gAz, we can define the 
Fourier transform as {z k) 



'qAz{Az) dAz, 



(44) 



which is often called the characteristic function of q^z- Considering then the 
expression 
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k — J k — 

= I Az''qAzi^z)dAz, (45) 

for n = 0, 1, 2, 3..., we see that, because q^z is a pdf, the last expression is the 
expectation value (A^") of A^", the so-called nth moment, and we have 

^"5fe"gA.(0) = (A^"). (46) 



In particular, it always holds that (A^;") = 1, since q^z is a pdf that is 
normalized to one {(Az^) = J q^ziAz) dAz = 1). Furthermore, (A^-*^) = 
/ AzqAz(yAz) dAz is the mean value of qAz- 

In the use of Fourier and Laplace transforms for solving the CTRW equations, 
we will concentrate on the asymptotic, large \z\ regime, which corresponds to 
small values of k (small wave-numbers correspond to large length-scales or 
wave-lengths). We thus can make a Taylor expansion of qAz{k) around k = 
and keep only a few low order terms, 

QAzik) = qAziO) + dkqAziO)k + ^dlqAziO)^ + (47) 



With Eq. fl46|) . the derivatives can be replaced with the moments, 

qAz{k) = l-i{Az)k-^{Az^)e + ... (4J 



(with (A^;*^) = 1). The Taylor expansion of a pdf in terms of the moments is 
practical because the moments are the natural characteristics of a pdf. Often, 
the distribution of spatial increments is assumed to be symmetric around 
z = 0, so that (Az) = 0, and the Taylor expansion writes in this case as 

qAz{k) = l-l{Az')e + ... (49) 



Temporal distributions, such as the time-step distribution qAt{At) in the wait- 
ing model, have the characteristic to be 'one-sided', i.e. they are defined and 
used only for t > 0, so that it is more appropriate to use Laplace transforms 
in this case, defined as 

oo 

qAtis) = J e-''^'qAt{At)dAt.. (50) 
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Straightforward calculations and the use of the respective definitions leads to 
the analogue of Eq. (146|) for Laplace transforms, 



oo 

-dsTQAtis) ^ = I Arq^tiAt) dAt = (At"), 

s=0 J 



(51) 



where the (At^) are again the moments. As with respect to z, we will focus 
on the asymptotic, large t regime, which corresponds to small values of s, and 
we make a Taylor-expansion of the Laplace transform around s = 0, replacing 
through Eq. (!5T!) the derivatives by the moments, 

qAt{s) = l-{At)s + ..., (52) 

in complete analogy to Eq. (1481) ((At°) = 1 is the normalization of q^t)- The 
Taylor expansions in Eqs. (HSj) . dlSD, and (15^ can of course only be used if 
the involved moments are finite. 



5.3.3 The symmetric and the one-sided Levy distributions 
The symmetric Levy distributions are defined in Fourier space as 

QazW = exp(-a|fc|"), (53) 



with < a < 2. It is not possible to express them in closed form in real 
space, with two exceptions, the case a = 2 is the usual Gaussian distribution 
(the Fourier back-transform of a Gaussian is a G aussian) , and t he case a = 
1 is known as the Cauchy distribution (see e.g. iHughed . Il995l . Chap. 4.3). 



As mentioned in Sec. 15.11 the Levy distributions for a < 2 all have infinite 
variance, and for a < 1 they even have an infinite mean value, so that the 
expansion in the form of Eq. (H^ is not applicable. We can though directly 
expand the exponential in Eq. (13^ and find the small k expansion as 



q^A:{k) = l-a\kr + 



(54) 



The case a = 2 corresponds obviously to the classical case of Eq. fj49l) with 
finite second moment (a = (1/2)(A2;^)) and zero mean (we are using only the 
symmetric Levy-distributions). 

If one-sided distributions with infinite variance are needed, then they can also 
be defined via Fourier space as a specific asymmetric Levy distribution, or. 
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more convenient for our purposes, as shown in iHughed (119951 ) (Chap. 4.3.2), 
they can be defined in Laplace space as 



exp(— &s^). 



(55) 



with b strictly positive, and where now < (3 < 1. These one-sided Levy 
distributions decay as t~^~^ for t oo, and they have the small s expansion 

qLhP(s) = l~bs^ + ... (56) 



Note that for /? = 1 we recover the finite mean case of Eq. (1521) . with b = (At), 
and the back-transform of Eq. (|55|) in this case yields the distribution 6{At—b), 
i.e. the time-steps are constant and equal to b {At = (At) = b). 

In the following, we will use the (small k) Fourier expansion for qAz{k) in the 
form, 

qAz{k) = l-a\k\'' + ..., (57) 



which for a < 2 corresponds to the Levy case, Eq. (!54|) . and for a = 2 it 
recovers the normal, finite variance case of Eq. (jUj), with a = {l/2){Az^). 
Correspondingly, the (small s) Laplace expansion for qAt{s) will be used in 
the form 

qAt{s) = l-bs^ + ..., (58) 



which for /5 < 1 yields the Levy distribution, Eq. (ISBjl . and for /? = 1 the 
normal, finite mean case of Eq. (1521) . with b = (At). 



5.3.4 Solving the CTRW equations with Fourier and Laplace transforms 

In this section, we will illustrate the use of the Fourier and Laplace (F-L) 
transform to solve the CTRW equation on the example of the waiting model, 
Eqs. (HUj) and (HT]) . The use of the respective convolution theorems, namely 

lf{x-y)g{y)dy-^f{k)g{k) (59) 

for Fourier transforms, and 

t 

J (pit - T)ip{T) dr ^ 4){s) ipis) (60) 
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for Laplace transforms (with /, g, 0, and ip any transformable functions), al- 
lows in the case of the waiting model to determine the solution in Fourier 
Laplace space: For the initial condition P{z,t = 0) = 6{t)6{z) and in the ab- 
sence of any source {S = 0), Eq. (l40l) turns into Q{k, s) = Q{k, s) qAz{k) qAt{s) + 
1, and Eq. fHT]) takes the form Pw{k,s) = Q{k, s)^w{s), and we can elim- 
inate Q and solve for P. Noting further that ^w{t) = Jt^ QAti^t) dAt = 
1 — /4 9At(At) dAt, so that ^wis) = (1 — qAtis))/s, we find 

1 - qAt{s) 



Pw{k,s) 



s[l- qAz{k)qAt{s)]' 



(61) 



which is known as the Montroll- Weiss equation (e.g. 


Montroll & Weiss. 


1965 


Zumofen & Klafter. 


1993; 


Klafter. Blumen & Shlesingei. 


1987 


)• 



Looking for asymptotic solutions, we insert the general form of the transformed 
temporal and spatial step distribution, Eqs. fl571) and fl58l) . respectively, into 
Eq. fl6T|) . which yields 



Pwik,s] 



bs^ + a\k\ 



(62) 



Unfortunately, it is not possible to Fourier and Laplace back-transform Pw{k, s) 
analytically. We can though use Eq. fH6l) for n = 2, namely 



{z\s)) = -dlPw{k = 0,s) 



(63) 



to determine the mean square displacement in the asymptotic regime (note 
that we set = at the end, which clearly is in the large \z\ regime). Inserting 

Pw{k, s) from Eq. (1^ into Eq. (jMD, without yet setting k = 0, we find 



is)) 



2^,21 7,1 2a-2 



2a'a'\k\ 



^^2/3+1 



+ 



aaia 



l)\k\ 



(64) 



The first term on the right side diverges for a < 1, and the second term di- 
verges for a < 2, so that {z^{s)) is infinite in these cases. This divergence 
must be interpreted in the sense that the diffusion process is very efficient, 
so that in the asymptotic regime, Pw has already developed so fat wings at 
large \z\ (power-law tails) that the variance, and with that the mean square 
displacement , of Pw is infinite, Pw has already become a Levy ty pe distribu- 
tion (see also lKlafter. Blumen fc Shlesingerl . 119871 : lBalescul . l2007al ). Of course, 
with the formalism we apply we cannot say anything about the transient 
phase, before the asymptotic regime is reached. We just note here that the 
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velocity model (Eq. (H2l) )m this regard is not so over-efficient, it allows super- 
diffusion with a more gradual build-u p of the fat wings of the distribution Py 
(iKlafter. Blumen fc Shlesingen . Il987l ). 



Less efficient diffusion can only be achieved in the frame of the waiting model 
for a = 2, i.e. for normal, Gaussian distributed spatial steps (see Eq. (1571) ). In 
this case, Eq. fl64|) takes the form 



is)) 



(65) 



(a = (1/2) (A^^) for a = 2). This expression is valid for small s, and with 
the help of the Tauberian theorems, which relate the power-law scaling of a 
Lap lace transforrn at snaall s t o the scaling in original space for large t (see 
e.g. iHughesl . Il995l : IPellei] . Il97lh it follows that 



{Z'is)))r^t 



/3 



(66) 



With our restriction < /3 < 1, diffusion is always of sub-diffusive character, 
and for /3 = 1 it is normal, as expected, since we have in this case waiting 
times with finite mean and variance (see Eq. flSSjl ). 



5.4 Including Velocity Space Dynamics 



Above all in applications to turbulent systems, and mainly to turbulent or 
driven plasma systems, it may not be enough to monitor the position and the 
timing of a particle, since its velocity may drastically change, e.g. if it interacts 
with a local electric field generated by turbulence. An interesting extension of 
the standard CTRW for these cases is to include, besides the position space 
and temporal dynamics, also the velocity space dynamics, which allows to 
study anomalous diffusive behaviour also in energy space. 

To formally define the extended CTRW that also includes momentum space, 
we keep Eq. (|9]) and Eq. (l39l) for position and time evolution as they are, and 
newly the momentum (or velocity) also becomes a random, dynamic variable, 
with temporal evolution of the form 

Pn = Apn + Apn-l + Apn~2 + ••• + ^Pl + Po, (67) 



with po the initial momentum, and the Api the momentum increments. Again, 
one has to specify a functional form for the distribution of momentum incre- 
ments qAp{Ap) in order to specify the random walk problem completely. The 
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solution of the extended CTRW is in the form of the distribution P{f,p, t) for 
a particle at time t to be at position r and to have momentum p. 



T he extended CTRW can b e tre ated b y Mont e - Carl o simulations, as done 



m 



i (2004 


), or in 


Islikei 


( 


2008) 



the extended CTRW has been introduced, which basically is a generalization 
of Eq. fHOl) and Eq. fH2l) . and a way to solve the equations numerically is 
presented. 



6 Prom random walk to fractional diffusion equations 



The purpose of this section is to show how fractional diffusion equations nat- 
urally arise in the context of random walk models. The starting point here 
are the CTRW equations for the waiting model, Eqs. (HUj) and (HTj) . which in 
Fourier Laplace space take the form of Eq. (I^Tj) . and on inserting the small 
k and small s expansion of the step-size and waiting-time distributions, Eqs. 
fl571) and fl58l) . respectively, the waiting CTRW equation takes the form of Eq. 
fl62l) . with a < 2 and (3 < 1. Multiplying Eq. fl62|) by the numerator on its 
right side, 

Pw{k,s){bs'^ + a\k\'') = bs^-\ (68) 



and rearranging, we can bring the equation for Pyy to the form 

a 
b 



s^Pwik, s) - = —Ikl^Pwik, s). (69) 



It is illustrative to first consider the case of normal diffusion, with /3 = 1 and 
a = 2, where according to Eqs. fl57|) and fl58l) we have a = (1/2) (A^^) and 
b = (At), so that 

s~Pw{k, s)-s' = J^^lkfp^ik, s) (70) 



Now recall how a first order temporal derivative is expressed in Laplace space, 
A^(^) ^ si^is) - s'm, (71) 

and a how a spatial derivative translates to Fourier space. 



dz' 



:f{z) {-ikrf{k) (72) 
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Obviously, for Pw{z,t = 0) = S{z), Eq. (170|) can be back-transformed as 



d,PUz,t) = ^^d^Pwiz^t), (73) 



so that we just recover the simple diffusion equation of the normal diffusive 
case. 



6.1 Fractional derivatives 



Fractional derivatives are a generalization of the usual derivatives of nth order 
to general non-integer orders. There exist several definitions, and in original 
space {z or t) they are a combination of usual derivatives of integer order 
and integrals over space (or time). The latter property makes them non-local 
operators, so that fractional differential equations are non-local equations, as 
are the CTRW integral equations. For the following, we need to define the 
Riemann-Liouville left-fractional derivative of order a, 



with r the usual Gamma-function, a a constant, n an integer such that n — l< 
a < n, a a positive real number, and / any suitable function. Correspondingly, 
the Riemann-Liouville right-fractional derivative of order a is defined as 

fiz) - J-D:^^ [ /(^O dz' (75) 



It is useful to combine these two asymmetric definitions into a a new, sym- 
metric fractional derivative, the so-called Riesz fractional derivative, 

' 2cos(7ra/2) 



The Riesz fractional derivative has the interesting property that its represen- 
tation in Fourier space is 

^^^Dy{z)^-\krf{k). (77) 



Comparison of this simple expression with Eq. (!72l) makes obvious that the 
Riesz derivative is a natural generalization of the usual derivative with now 
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non- integer a. 

To treat time, a different variant of fractional derivative is useful, the Caputo 
fractional derivative of order /3, 

^^^^D^Mt) = — / , ^(t') dt', (78) 



with n an integer such that n — 1 < [3 < n, and any appropriate function. 
The Caputo derivative translates to Laplace space as 

(^)D^^{t) ^ s^tPis) - s^-V(O), (79) 



for < f3 < 1, which is again a natural generalization of Eq. flTTl) for the usual 
derivatives for now non-integer j3 (the Caputo derivative is also defined for 
/? > 1, with Eq. fl7^ taking a more general form). 



Furth er details about fractional deriv atives can be fo und e.g. in iPodlubny 



(Il999l ) or in the extended Appendix of iBalescul (l2007bl ). 



6.2 Fractional diffusion equation 



Turning now back to Eq. fl69l) . we obviously can identify the fractional Riesz 
and Caputo derivatives in their simple Fourier and Laplace transformed form, 
Eqs. f l77|) and fl7^ . respectively, and write 

^""^D^Pwiz, t) = ^ ^""^D^^^P^iz, t) (80) 



/^From this derivation it is clear that the order of the fractional derivatives, a 
and P, are determined by the index of the step-size (qaz) and the waiting time 
{lAt) Levy distributions, respectively. It is also clear that Eq. (I5UI) is just an 
alternative way of writing Eq. (1^ or (1^ . and as such it is the asymptotic, 
large \z\, large t version of the CTRW equations fHOj) and fHT]) . It allows though 
to apply different mathematical tools for its analysis that have been developed 
specially for fractional differential equations. 

As an example, we may consider the case /? = 1 and < a < 2, where the 
diffusion equation is fractional just in the spatial part, 

dtPUz,t) = ^ (''^D^,^Pw{z,t) (81) 
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In Fourier Laplace space, this equation takes the form 

which, on applying the inverse Laplace transform, yields 

Pw{k,t) = exp(^-^\k\''ty (83) 

which is the Fourier-transform of a symmetric Levy-distribution with time as 
a parameter (see Eq. (155]) ). and with index a equal to the one of the spatial 
step distribution g^z- Thus, for a < 2, the solution has power-law tails, and 
the mean square displacement (or variance or second moment) is infinite, as 
we had found it in Eq. flM|) . For a = 2, the solution Pw{z,t) is a Gaussian 
(the Fourier back-transform of a Gaussian is a Gaussian), and we have normal 
diffusion. 



7 Action diffusion in Hamiltonian systems 



So-far, our starting point for modeling diffusion was mostly the random walk 
approach and a probabilistic equation of the Chapman Kolmogorov type (Eq. 
( pTl) ). Here now, we turn to Hamiltonian systems, and we will show how from 
Hamilton's equations a quasi-linear diffusion equation can be derived. This 
diffusion equation is of practical interest when the Hamiltonian system consists 
in a large number of particles, so that it becomes technically difficult to follow 
the individual evolution of all the particles. 

Let us consider a generic iV— degrees of freedom Hamiltonian system with 
Hamiltonian H (q, p) and equations of motion given by 



It 
dp 
Itt 



dH 
dp ' 
dH 

9q ' 



B4) 
S5) 



where q = (gi, ...,qN) and p = (pi, ...,Pn) are the canonical coordinates and 
momenta, respectively. In order to be integrable such a system should have 
A'^ independent invariants of the motion, co rresponding to an equal number of 
symmetries of the system (Goldstein, 19801 ). The integrability of a system is a 
very strong condition which does not hold for most systems of physical interest. 
However, in most cases we can consider our system as a perturbation of an 
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integrable one and split the Hamiltonian accordingly to an integrable part and 
a perturbation. Then the description of the system can be given in terms of the 
action-angle variables of the integrable part (Note that a periodic in t egrab le 
system can always be transformed to action- angle variables (lGoldsteinl . ll980l )). 
so that we can write 



H{J,e,t) = Ho{3) + eH^{3,e,t), 



with J = (Ji, J/v) and = (6*1, ...9^) being the action and angle variables, 
respectively. Hq is the integrable part of the original Hamiltonian and Hi is 
the perturbation. The parameter e is dimensionless and will be used only for 
bookkeeping purposes in the perturbation theory; it can be set equal to unity, 
in the final results. The evolution of the integrable system Ho is given by the 
following equations of motion 



J = 

e=u}ot + Oo, 



^7) 



where cjq = dHo/dJ are the frequencies of the integrable system Hq. The N 
action variables correspond to the N invariants of the motion required for the 
integrability of the system. 

The perturbation Hi leads to the breaking of this invariance due to its 
dependence. The derivation of a quasilinear diffusion equation in the action 
space is the subject of this section, and the method to be used is the canonical 
perturbation theory applied for finite time intervals. This method of deriva- 
tion is based on first principles and does not imply any statistical assumptions 
for the dynamics of the system, such as the presence of strong chaos resulting 
in phase mixing or loss of memory for the system. Moreover, it is as sys- 
tematic as the underlying perturbation scheme of the canonical perturbation 
theory, so it can be extend ed to higher or der and provide results beyond the 
quasilinear approximation (jKominisi . 120081 ). Also, it is important to note that 
the method makes quite clear what physical effects are taken into account in 
the quasilinear limit and what effects are actually omitted. It is worth men- 
tioning that non-quasilinear diffusion has been studied both analytical l y and 



numerically for a variety of physical systems (Carv. Escande &: Vergal . 



Helander fc Kiellberd . Il994 
1999h . 



Benisti fc Escande. 1998: Laval fc Pesme 



1990 



1983 



The basic idea of the canonical perturbation theory is the search of canonical 
transformations for the perturbed system (i.e. transformations which preserve 
the Hamiltonian structure of the system) under which the new (transformed) 
Hamiltonian is a function of the action only. For a near-integrable system this 
can be done approximately, and the new actions correspond to approximate 
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invariants of the motion which contain all the essential features of the phase 
space structure. The transformations involved in canonical perturbation the- 
ory are expressed in terms of the so-called mixed- variable generating functions. 
These can be fun ctions of a subse t of the old variables along with a subset 
of the new ones (iGoldsteinl . Il980l ). Thus, the transformation from (J, ^) to 
(J, ^) can be expressed by a generating function of the form S{J,6,t). The 
transformation equations are given in the following implicit form: 



J = J + e 

e = e + e 



dS{3,0,t) 

80 

dS{3,e,t) 

dJ ■ 



(89) 
(90) 



Following a standard procedure ( Goldstein . 19801 ). we seek a transformation 
to new variables (J, 6) for which the new Hamiltonian H is a function of the 
action J alone. Expanding S and H in power series of a small parameter e 



(91) 
(92) 



where the lowest-order term has been chosen to generate the identity trans- 
formation J = J and = 0. The old action and angle can be also expressed 
as power series in e: 



" — J~r6 \ € ~r 



= 



do dO 

dSi{3,0,t) ^dS2{3,0,t) 



dJ 



+ ..., 



and the new Hamiltonian is 



H{J,0,t) = H{3,0,t) + 



03(3,0, t) 

di ■ 



(93) 
(94) 

(95) 



By substituting the respective power series in Eq. (|95|) and equating like powers 
of e for the zero order we have 



Hq = Ho, 



(96) 



while in the first and second order we have the equations 



'dt °~dO ^' 



1,2, 



(97) 
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with 




+ 



dHi dSi 
dJ dO 



(98) 
(99) 



providing the first and second order generating function 5*1 and 5*2, respec- 
tively. The latter are linear partial equations which can be solved in a time 
interval of interest [to, t] by the method of characteristics (i.e. integration along 
the unperturbed orbits). Note that Hi and H2 are arbitrary functions which 
can be set equal to zero (for the application of the canonical perturbation the- 
ory in infinite time i ntervals, these fu nctions have to be chosen so that they 
cancel secular terms f Goldstein . 19801 )). For a general perturbation of the form 



Hi 



EH. 



,6 



(100) 



the solution for the first order generating function 5*1 can be written as 

t 

Si{J, 6, t; to) = - E e^'"^^-^"*) / H^J, s)e^-^o^(is. (101) 



to 



Similarly, the solution for S2 can be readily obtained. The resulting expression 
(too lengthy to be presented here) is a periodic function of 0. This is the only 
information we need for 5*2, since its exact form will not be involved in our 
calculations. 



In the following, we show that the results of first order perturbation theory can 
be utilized in order to provide an evolution equation for the angle-averaged 
distribution function, which is accurate up to second order with respect to 
the perturbation parameter e, namely a quasilinear action diffusion equation. 
Therefore, we can relate results from perturbation theory applied for a sin- 
gle particle motion, to the distribution function, describing collective particle 
motion. The latter is of physical interest in all cases where a large number of 
particles is involved in collective phenomena so that a statistical approach is 
required. 



The evolution o f the phase space distribution function F is governed by Liou- 
ville's equation (IGoldsteiru . Il980l ) 



dF 



0, 



(102) 



32 



where [., .] denotes the Poisson bracket, defined as [/i, = Vq/i-Vp/2— Vp/i- 
Vq/2 with q and p being the canonical positions and momenta, respectively. 
This equation simply expresses the incompressibility of the Hamiltonian fiow 
and the invariance of the number of particles. It is well-known that for an 
integrable system, any function of the invariants of the motion (actions) is a 
solution of the Liouville's equation. For the case of a near-integrable system, 
an approximate distribution function can be obtained as a function of the 
approximate invariant of the motion, namely the new actions J, so that we 
can write 



F{3,e,t) = F{3), 



(103) 



with J given implicitly by Eq. flU51) . To second order, with respect to e we have 

e2 d 



3 = 3- eAiJ + ■ (AiJAiJ) + e'AsJ, 



(104) 



with 



AiJ(J,6>,t;to) 



dS,{3,e,t;to) 
dd 



1,2. 



(105) 



Substituting (11041) in (I103P and utilizing a Taylor expansion with respect to e 
we have 



F(J,.,t)^F(J)-e^.(A,J)4| 



(AiJAiJ) 



dF{3) 



83 



(106) 



Noting that AiJ(J, 6, to] to) = we have -F(J, 6, to) = F{3), and by averaging 
over the angles, we obtain 



f{3,t) = f{3,to) + 



2 03 



((AiJA.J))^ • 



(107) 



where / is the angle-averaged distribution function 
/(J,t) = (F(J,0,i))0, 



(108) 



and we have used the fact that (AjJ)^ = as obtained from Eqs. (11051) and 
the fact that Si are sinusoidal functions of 6 . Taking the limit t ^ to in Eq. 
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(11071) . we finally obtain the quasilinear action diffusion equation 



dt 



9J 



D(J,t) 



dl 



(109) 



with 



l g((AiJA,J))g 
^^^'*^"2 dt 



(110) 



being the corresponding quasilinear diffusion tensor. Identifying that Ai J cor- 
responds to the first order action variation, we see that D in Eq. (II 101) cor- 
responds to the common definition of the diffusion tensor, as provided by the 
statistical approach, which is base d on the Kra r ners-M oyal expansion of the 
master equation (see Sec. 13.41 and Ivan Kampenl . Il98ll ). 



It is worth mentioning that this form [Eq. (I109p ] of the action diffusion equa- 
tion is similar to the diffusion equation obtained with the utilization of the 
Fick's Law [Eq. (pTj) ]. However, it has been shown that for any Hamiltonian 
system, this form of the action diffusion equation is equivalent to the Fokker- 
Planck equation [Eq. due to the fact that the corresponding drift velocity 
and diffusion terms are related through 



;iiii 



where z = J for the case of th e action diffusion equation (see Chap. 4 in 



Lichtenberg fc Liebermanl . Il992l . and references therein) . This property is 
implied directly from the canonical form of the underlying equations of motion. 
More generally, it has been shown that this form of the diffusion equation has 
further relat ion with the more g e neral class of the microscopically reversible 



systems (see iMolvig fc Hizanidid . Il984j . and references therein). 



It is important to emphasize that the derivation procedure described here does 
not prerequisite any statistical assumption and is only based on the under- 
lying equations of motion. In contrast to other statistical approaches, it is 
not necessary to assume strong stochasticity related to the completely chaotic 
regime where loss of memory takes place and the orbits are completely decor- 
related. Therefore, Eq. (I109p . is capable of describing not only diffusion in the 
almost homogeneous "chaotic sea" in the phase space, but also intermittent 
motion in an inhomogeneous phase space structure where resonant islands 
and chaotic areas are interlaced. It is worth mentioning that the diffusion ten- 
sor D is time-dependent. The dependence of D on the actions, for the case 
of time-periodic perturbations, is through smooth l ocalized f unctions which 
are centered around the corresponding resonances (IKominisl . [2OO8,). In the 
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limit t ^ oo these localized functions tend to Dirac delta functions, which 
commonly appear in standard derivations of the quasilinear diffusion tensors 



( IZaslavski fc Filonenkd . Il968l ). However, the consideration of this limit in the 
derivation of D corresponds to an extension to infinity of the limits of in- 
tegration in the derivation of 5*1, in Eq. fllOll) . which implies an ergodicity 
assumption and a steady-state approach. In the most general case, the diffu- 
sion tensor D [Eq. (11 101) ] is capable of describing not only steady-state but also 
transient diffusion phenomena, through its time dependence. In this case the 
time scales of the action diffusion are determined by the time-dependence of D 
as well as by the factor in the r.h.s. of Eq. (11091) . which implies an actually 
slow diffusion process. Note that these time scales come naturally into play 
and there is no need for a priori separation of the distribution function in slow 
and fast varying parts as in many heuristic derivations of the Fokker-Planck 
equation. 

In the previous paragraphs we have derived a general quasilinear action diffu- 
sion equation for a Hamiltonian system, with a minimum of assumptions on 
the underlying dynamics. The fact that this method of derivation is closely 
related to a systematic and rigorous perturbation scheme allows for extending 
these results beyond the quasilinear approximation. Therefore we can carry 
our perturbation scheme to higher order in order to provide a hierarchy of 
diffusion equations having higher-order derivatives of the action distribution 
function with respect to the action, in direct analogy to statistical approaches 
whe re higher-ord e r Kra mers- Moyal expansions are considered (see Sec. 13. 4[ 



and Ivan Kampenl . Il98ll ). Note that the most appropriate method for handling 



calculations inv olved in higher -order perturbation theory is the utilization of 



Lie transforms (jKominisi . l2008l ) . The hierarchy of higher-order diffusion equa- 
tions does not only provide better accuracy with respect to the perturbation 
parameter e, but is also capable of describing non-Gaussian evolution of the 
distribution function and resonant processes between the particle and beats 
of multiple spectral components of the perturbation, known as nonlinear res- 
onances. 



8 Other ways to model anomalous diffusion 



Escande fc Sattin fl2007h review and discuss under what circumstances the 
Fokker-Planck equation (Eq. (1341) ) is able to model anomalous diffusion. In 
summary, the FP equation, which is a local model, is able to model anomalous 
diffusive behaviour in cases where there is a non-zero drift velocity, V{z) ^ 0, 
anomalous diffusion is thus based on drift effects. 



Klafter. Blumen fc Shlesinged (119871 ) briefly review attempts of using the Fokker- 
Planck equation with zero drift velocity, but spatially or temporally dependent 
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diffusion coefficient. In order to account for anomalous diffusion, the spatial 
or temporal dependence of the diffusion coefficient must though be chosen in 
very particular ways, which are difficult to interpret physically. 



Lenzi. Mendes &: TsallisI (120031 ) shortly discuss non-linear diffusion equations, 
which have the form of the simple diffusion equation as in Eq. (130!) . with P 
though raised on one side of the equation to some power 7. 



9 Applications in Physics and Astrophysics 



CTRW has successfully been applied to model various phenomena of anoma- 
lous diffusion, including sub- and super-diffusive phenomena, in the fields of 
physics, chemistry , astronomy , biology, and economics (see the references in 
Metzler fc Klafterl . l2000l . l2004j i . 



Laboratory plasma in fusion devices (tokamaks) show a variety of anomalous 
diffusion phenomena. iBalescu fll995h was the first to apply CTRW to plasma 



physic al problems. Later. Ivan Milligen. Sanchez fc Carrerad ( 120041 ) and lvan Milligen. Carreras &: Sanch 
( I2OO4J ) developed a CTRW model for confined plasma, the critical gradient 
model, which was able to explain observed anomalous diffusion phenomena 
such as 'up - hiir t ransport, where particles diffuse against the driving gradi- 
ent. Ilslikerl 



Islike i (120081) studied the same physical system, with the use though of 



the extended CTRW that includes momentum space dynamics, and they stud- 
ied the evolution of the density and temperature distribution and the particle 
and heat diffusivities. 



Also in astrophysical plasmas anomalous diffusion is ubiquitous, there are 
many astrophysical systems where non-thermal (i.e. not Maxwellian distributed) 
particles are directly or indirectly, through their emission, observed. 



Dmitruk et al.l (l2003l . 120041 ) analyzed the acceleration of particles inside 3-D 
MHD turbulence. The compressible MHD equations were solved numerically. 
In these simulations, the decay of large amplitude waves was studied. After 
a very short time (a few Alfven times), a fully turbulent state with a broad 
range of scales has been developed (Fig. 15]). 

The magnetic field is directly obtained from the numerical solution of the 
MHD equations, with electric field derived from Ohm's law. It is obvious that 
the electric field is an intermittent quantity with the high values distributed 
in a less space filling way. Magnetic and electric fields show a broad range of 
scales and high degree of complexity. The energy spectrum of the MHD fields 
is consistent with a Kolmogorov-5/3 power law. The structure of the velocity 
field and the current density along the external magnetic field (J^) can be 
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Fig. 5. Visualization of the turbulent magnetic field | B \ (top) and electric field | E \ 
(bottom) in the simulation box. High values are in yellow (light) and low values in 
blue (dark). 

seen in Fig. [6l The formation of strong anisotropies in the magnetic field, the 
fluid velocity and the associated electric field is observed. The overall picture 
is that current sheet structures along the DC field are formed as a natural 
evolution of the MHD fields. 




Fig. 6. Cross section of the current density Jz along the external magnetic field in 
color tones. Yellow (light) is positive Jz, blue (dark) is negative, and the superposed 
arrows represent the velocity field. 

Following thousands of particles particles inside the simulation box, we can 
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learn many of the statistical properties of their evolution, e.g. the mean square 
displacements \J < Ax^ >, a/ < Av"^ >, or the velocity distribution etc. can be 
determined. Electrons and ions are accelerated rapidly at the nonlinear small 
scale structures formed inside the turbulent volume, and non-thermal tails of 
power-law shape are formed in the velocity distributions. Most particles seem 
to escape the volume by crossing only a few of the randomly appearing current 
sheets. A few particles are trapped in these structures and accelerated to 
very high energies. The Fokker-Planck equation is not the appropriate tool to 
capture particle motion in the presence of the random appearance of coherent 
structures inside such a turbulent environment. 



Vlahos. Isliker fc Lepretil (120041 ) performed a Monte Carlo simulation of the 
extended CTRW in position and momentum space, in application to flares in 
the solar corona, with particular interest in the appearance of the non-thermal 
energy distributions of the so-called solar energetic particles. 



10 Summary and Discussion 



Brownian motion is a prototype of normal diffusion, and its analysis has 
brought forth a number of tools that today are very much in use for modeling 
a wide variety of phenomena. Normal diffusion occurs in systems which are 
close to equilibrium, like the water in Brown's experiment. It has now become 
evident that phenomena of anomalous diffusion are very frequent, because 
many systems of interest are far from equilibrium, such as turbulent systems, 
or because the space accessible to the diffusing particles has a strange, e.g. 
fractal structure. The tools to model these phenomena, continuous time ran- 
dom walk, stochastic differential equations, and fractional diffusion equations, 
are still active research topics. 
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